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Abstract 

H2 is the most abundant interstellar species. Its deuterated forms (HD and D2) are also 
significantly abundant. Huge abundances of these molecules could be explained by 
considering the chemistry occurring on the interstellar dust. Because of its simplicity, 
Rate equation method is widely used to study the formation of grain-surface species. 
However, since recombination efficiency of formation of any surface species are heav¬ 
ily dependent on various physical and chemical parameters, Monte Carlo method would 
be best method suited to take care of randomness of the processes. We perform Monte 
Carlo simulation to study the formation of H2, HD and D2 on interstellar ices. Ad¬ 
sorption energies of surface species are the key inputs for the formation of any species 
on interstellar dusts but binding energies of deuterated species are yet to known with 
certainty. A zero point energy correction exists between hydrogenated and deuterated 
species which should be considered while modeling the chemistry on the interstellar 
dusts. Following some earlier studies, we consider various sets of adsorption energies 
to study the formation of these species in diverse physical circumstances. As expected, 
noticeable difference in these two approaches (Rate equation method and Monte Carlo 
method) is observed for production of these simple molecules on interstellar ices. We 
introduce two factors, namely, S / and [i to explain these discrepancies: 5 / is a scal¬ 
ing factor, which could be used to correlate discrepancies between Rate equation and 
Monte Carlo methods, p factor indicates the formation efficiency under various cir¬ 
cumstances. Higher values of p indicates a lower production efficiency. We found that 
P increases with a decrease in rate of accretion from gas phase to grain phase. 
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1. Introduction 

Molecular hydrogen is the most abundant and simplest species in the Interstel¬ 
lar Medium (ISM). Indeed, H2 is the most important molecule, as it is the first and 
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foremost precursor to create other complex molecules. Simpler molecules are mainly 
formed on grain surfaces and then desorbed to gas phase (Gould & Salpeter, 1963). Im¬ 
portance of grain chemistry were already described in the earlier literatures (Hasegawa, 
Herbst & Leung, 1992; Chakrabarti et. al., 2006ab; Cuppen & Herbst, 2007; Das et al., 
2008ab; Das et al., 2010, 2013ab; Das & Chakrabarti, 2011; Majumdar et al., 2012; 
2013; Das et al. 2014). Sometimes even without explicit grain chemistry, effective 
rates were used in studying the formation of bio-molecules (Chakrabarti & Chakrabarti, 
2000ab) For large grains or high accretion rates. Rate equation method (Biham et al., 
2001) can be used to explain abundances of surface species. But due to inherent ran¬ 
domness of the processes, especially when grain size is small and/or accretion rate 
is small, simple Rate equation method is not good enough to explain recombination 
efficiency of interstellar surface species. It is then essential to consider Monte Carlo 
method to account all the features. 

Despite low elemental abundances of atomic deuterium (having D/H ratio of ~ 10“^ 
according to Linsky et al., 1995), several complex molecules are found to be heavily 
fractionated (Majumdar, Das & Chakrabarti, 2014ab) in ISM. As like the normal hy¬ 
drogen molecules, its deuterated forms could also be synthesized on interstellar grains. 
Chakrabarti et al. (2006ab) carried out Monte Carlo simulation to find out exact re¬ 
combination efficiency for the formation of H 2 . Following similar approach, here we 
use both Rate equation method and Monte Carlo method to find out the various aspects 
during the formation of H 2 , HD and D 2 . 

The plan of this paper is the following. In Section 2, computational details are 
presented. Implications of results are discussed in Section 3. Finally, in Section 4, we 
draw our conclusions. 


2. Computational details 


2.7. Rate Equation Method: 

A small chemical network containing H and D is used here. Let at any time f, Nx 
be the number of species x on a grain having S number of adsorption sites. Governing 
equations for all surface species in our network could then be written as follows: 


dNn 

dt 


- Fff - WhNh - 2auN]jlS - aHoNuNolS, 


dNn, 

dt 


- Fhi + RHiOnN^j is - WH2NH2, 


dND2 

dt 


— Fd 2 + RD2t^DN^IS - B/jjiYDj, 


( 1 ) 

( 2 ) 

( 3 ) 


dNo 

dt 


- Fo — WdNd — ^aoN^jS - aHoNMNolS, 


( 4 ) 


dN MD 

— -7— - Fhd + RhdOhdNhNdIS - WhdNhd, ( 5 ) 

dt 
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where, F^, Ox and Wx represents accretion, hopping and desorption rates respectively 
for species x in the units of sec~^. jix represents spontaneous desorption factor. It is 
assumed that (1 - /i^) factor of species could go to gas phase just after its formation 
on grain. From experimental findings of Katz et al. (1999), we use = 0.33 for 
olivine grains and = 0.413 for amorphous carbon grains. Here, we use similar 
spontaneous desorption factor for FID and D 2 as well. Accretion rate (Fj^) of x‘^ species 
is calculated from, 

Fx^txA{V)Nx, ( 6 ) 

where, is the sticking coefficient, A is the surface area of the grain in the units of crn^, 
< y > is average thermal velocity (< V >= cm sec”') and rix is the gas phase 

concentration of x'^ species (in cmT^). 


2.1.1. Sticking coefficient 

Since the sticking coefficient of //2 is ~ 0 at F = 10 K (Leith & Williams, 1985), 
here. We consider, sticking coefficients (tx) of H 2 , HD and D 2 to be ~ 0. Thus, from 
Eqn. 6, we have FFfjD & Fd, = 0. In reality, sticking probability of a species is 
mainly governed by kinetic energy of the incoming species and adsorption energy of 
that species with grain surface. Since, hydrogen is the lighter species, hydrogenation 
reaction is the fastest reaction on the grain surface. In order to determine the sticking 
coefficient, various attempts were made over the last few decades. Buch & Zhang 
(1991) performed molecular dynamical simulation to numerically evaluate sticking of 
hydrogen atoms with cluster of water molecules. According to their study, sticking 
coefficient depends on the following: 

5 - (KbT/Eq + l)-2, (7) 


where, EqIKb — 102 K for H atom and 200 K for D atom. As per Chaabouni et al. 
(2012) sticking parameter of yflO/T could be adopted. Recently, Matar et al. (2010) 
performed an experiment to model the sticking parameters of H 2 and D 2 and its depen¬ 
dence on impinging molecular beam temperature. They found out an analytic formula 
(Eqn. 8) to describe sticking coefficient at any temperature T. Erom the outcome of 
their experiments and the interpretation of Chaabouni et al. (2012), sticking coefficient 
could be parameterized for hydrogen and deuterium on silicate surfaces under interstel¬ 
lar conditions. They pointed out that sticking coefficients of these species with silicate 
grains behave in the same as they would on icy dust grains. Their prescribed variation 
of sticking parameter is as follows: 


S(T)^So 


(1 +/3T/T0) 
d+T/Tor’ 


( 8 ) 


where, 5o is sticking coefficient of particles at zero temperature and Fq is a critical 
temperature obtainable using experimental data, cr defines the geometry of incident 
beam. Here, in our simulation, unless otherwise stated we always consider the sticking 
coefficient of all the species 1. Eor a special case, we consider the assumption of 
Chaabouni et al. (2012). 
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2.1.2. Accretion, Diffusion and Desorption 

An effective accretion rate on a typical grain is defined as, (p^ - 
where, '^ifgr^i is fraction of grain that could be occupied by x, number of species. 
Hopping rate of a species x is calculated by, 

= (9) 


where, Vx is vibrational frequency and Ehx is energy barrier for diffusion. Vibrational 
frequency of any species x is calculated by. 


Vx 


2sEdx 
n^nix ’ 


( 10 ) 


where, nix is mass of species x, Ed is adsorption energy of species x, s is surface density 
of sites in the unit of cnT^ (s = 10^"^ cm“^ is used). Finally, desorption rate of a species 
X is calculated by, 

Wx^Vx^M-EbxlKbT). ( 11 ) 


2.1.3. Binding energies 

Binding energies are important to decide the mobility of surface species. From 
experimental findings of Pirronello et al. (1997, 1999), Katz et al. (1999) concluded 
that atomic hydrogen moves much more slowly than what is normally used in vari¬ 
ous simulations. In our simulation, we use these experimental findings. In Katz et al. 
(1999), Eb and Ed were obtained for H atom. However, for H 2 molecules, only Ed was 
defined. Since interaction between adsorbed H and D atoms and grains are different, a 
zero point energy correction due to isotopic effects should to be considered because of 
which binding energies of deuterated species must differ. D atom being heavier than H 
atom, it would sit at a level lower than El in any potential surface if zero-point energy 
is included (Lipschitz, Biham & Herbst, 2004). Caselli et al. (2002) proposed a differ¬ 
ence of 2 meV between desorption energies and barriers against diffusion. Following 
this assumption, Lipshtat, Biham & Herbst (2004) (hereafter LBH) considered a 2 meV 
energy difference between H and D atoms for barriers against diffusion. They consid¬ 
ered energy barriers against diffusion for H atoms to be 35 meV whereas for D atoms 
they considered it is to be 37 meV. In case of desorption energies, they considered 
for H atoms it is 50 meV and for D atoms it is 60 meV. So, in case of barrier against 
desorption, LBH considered an energy difference of 10 meV. As per LBH, this larger 
difference in case of desorption was considered because barrier against diffusion bal¬ 
ances zero-point energy of a potential well with that of a saddle point while desorption 
does not possess a saddle point. 

In absence of a consensus on binding energies, we use three sets of energy barriers 
in our simulations. These sets are shown in Table 1. First set corresponds to experi¬ 
mental binding energies which were taken from Katz et al. (1999) for Olivine grains. 
Second set corresponds to experimental values obtained for Amorphous carbon grains 
(Katz et al., 1999) and third set is same as LBH. Binding energies for deuterated species 
are calculated by using following scaling relations. 


Eb(D) = Eb(H) X 


31 {Eb{D) from LBH) 
35 (Eb(H) from LBHY 


4 



Table 1: Binding energies used in our simulations 


Species 

Olivine grain 
(set 1) 

Amorphous carbon grain 
(set 2) 

Lipshtat, Biham & Herbst (2004) 
(set 3) 


Eh (in meV) 

Ed (in meV) 

Eh (in meV) 

Ed (in meV) 

Eh (in meV) 

Ej (in meV) 

H 

24.7 

32.1 

44.0 

56.7 

35.0 

50.0 

D 

26.1 

38.5 

46.5 

68.0 

37.0 

60.0 

H 2 

8.1 

27.1 

14.0 

46.7 

12.4 

41.2 

HD 

9.8 

32.5 

16.8 

56.0 

14.8 

49.4 

D 2 

9.8 

32.5 

16.8 

56.0 

14.8 

49.4 


Ed{D) = Ed{H) X 


60 {Ed{D) froniLBH) 
50 (Ed(H)fromLBHy 


Ed(HD,D2)^EdD)x 


46.7 (Ed(H 2 ) for amorphous carbon from Katz et al. (1999)) 

56.7 (Ed(H) for amorphous carbon from Katzet al. (1999)) 


Following Hasegawa, Herbst & Leung (1992); Das et al. (2008) and references therein, 
we consider Ed(x) = 0.3 Eh(x) for EI 2 , HD and D 2 molecules. In case of set 3, desorp¬ 
tion energy for H 2 was not available. We use following relation to compute desorption 
energy of H 2 for set 3: 


Ed{H2) = Ea{H) X 


46.7 (Ed(H 2 ) for amorphous carbon from Katzet al. (1999)) 

56.7 {Ed{H) for amorphous carbon from Katz et al. (1999)) 


In reality, nature of grain surface need not be simple. There might be various types 
of lattice defects and various phases. Defects with enhanced binding energies (Hol- 
lenbach & Salpeter, 1971) could allow formation of H 2 up to 50 K. In our simula¬ 
tion, we do not consider such types of grains and use a very low temperature window 
{%K-25 K). 


2.2. Procedures adopted in a Monte Carlo Simulation 

In order to preserve randomness in molecular formation process, it is essential to 
develop Monte Carlo algorithm to mimic exact scenario. Formation of molecules on 
grains are mainly due to two types of forces, namely, (i) the weak Van der Walls in¬ 
teraction and (ii) strong covalent bond (chemisorption). For low temperature region 
and/or for low surface temperature, atoms are weakly bound to surfaces via first type 
of force called physisorption. Here energy is of the order of a few meV. As species are 
weakly bound to surfaces, they can move around the surfaces via quantum mechanical 
tunneling or by thermal hopping. These processes are completely random. A species 
can move in any directions on the grain surface and could combine with another species 
to form a third species. We assume the surface to be a square lattice having nxn number 
of sites in each direction. To mimic spherical nature of the grain, we consider periodic 
boundary condition. Simulation on a complex shaped grain is outside the scope of the 
present work, and will be taken up elsewhere. 

Governing equations (Eqn. 1-5) are implemented in Monte Carlo method by as¬ 
suming a four step process. First step is accretion of gas phase species onto a grain 
surface. Second step is diffusion of surface species. Third step is reaction among hop¬ 
ping species (i.e., at chance-meeting site). Finally, fourth step is evaporation of the 
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product to gas phase. Reaction could occur via thermal diffusion or quantum mechani¬ 
cal tunneling. Since, we are considering experimental values of binding energies, only 
thermal hopping is considered. Production of surface species via Eley-Rideal mecha¬ 
nism (Das et al., 2008 and references therein) are also considered. All the possibilities 
in Monte Carlo method are handled by generating the random numbers as and when 
required. 

In our simulation, smallest time scale is the hopping time scale of H or D atoms. 
Accretion time scale is much longer than hopping time scale. Therefore, for simplicity, 
we are dropping a hydrogen atom after every QhIFh steps and a deuterium atom in 
every odIFo steps. Dropping locations are automatically chosen by generating a pair 
of random numbers. Depending on the number (n) of sites in each direction of the 
grain concerned, random numbers are scaled. While dropping FI or D, if the designated 
location is already occupied by H 2 , HD or D 2 then we look for a different site for its 
landing. However, if the site is occupied hy D or H then HDID 2 or H 2 IHD could be 
formed. Depending on binding energies, surface species are allowed to hopp towards 
any of the four directions. Directions are chosen by generating random numbers. When 
an atom meets another atom via thermal hopping, a new species could be formed with 
certainty (i.e., no activation barrier is considered). Evaporation of surface species are 
also handled by generating random numbers. Desorption rate for any surface species is 
calculated by using Eqn. 9 which implies that species x will leave after IjWx seconds. 
We are generating random numbers after every I/ah seconds to check this desorption 
probabilities. If this number is less than Wt/a^ then selected x species are released from 
grain phase. There are another type of desorption, namely spontaneous desorption, 
which are also handled by generating random numbers. This desorption term decides 
whether any newly formed species will remain on grain surface or desorbed from the 
surface. 


3. Results and Discussion 


In Eig. 1, chemical evolution of the surface species (H, H 2 , D, HD & D 2 ) are 
shown. Results from both Rate equation and Monte Carlo method are shown. Clear 
difference between these two methods are visible. Here, we consider amorphous car¬ 
bon grains having hydrogen (atomic form) number density {nu) - 10 ^ cnT^, T = \5 K 
and initial atomic DjH ratio (hereafter r^) = 0.1. Surface abundance of H atom has 
achieved a steady state after ~ 10^ seconds, whereas, for D atoms, steady state is 
achieved after ~ 10^ seconds. Eor this set of physical parameters, HD is the most 
abundant species on grain surface. A steady state would be achieved beyond 2x10^ 
seconds. It is to be noted for this choice of physical parameters HD and D2 are over¬ 
produced in the Rate equation method (Eig. 1). 

In Eig. 2, efficiency window for the formation of H 2 , D 2 and HD molecules are 
shown for Olivine grain (Eig. 2a), Amorphous carbon grain (Eig. 2b) and a type of 
grain having binding energies intermediate of Olivine and Amorphous Carbon type 
(Eig. 2c, with set 3 energy values of Table 1) (Eig. 2c) are shown for rin - 10^ 
and rj) - 0.1. Production efficiency is defined by following relation; 


t1h2 = 


2Rh2 
Fh ’ 


( 12 ) 
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Figure 1: Time evolution of abundances of H, H 2 , D, HD & D 2 by Rate equation method (dotted lines) and 
Monte Carlo method (solid lines). Simulation were carried out for an amorphous carbon grain kept at 15 K, 
rifi = 10^ cm“^ and = 0.1. 
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Efficiency 



Figure 2: Efficiency window for vaiious sets of binding energies for (a) Olivine, (b) Amoiphous Carbon and 
(c) for LBH case. Solid lines are results obtained from Monte Carlo method whereas dotted lines are results 
from Rate equation method. 
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(13) 


'702 = 


2Rd, 
Fd ’ 


tIhd 


Rhd 

Fh + Fd 


(14) 


where, Rd^ and Rhd are gas phase productions of H 2 , D 2 and HD respectively. 
Gas phase production of these species solely depends on barrier against desorption 
(i.e., thermal desorption). Thus, 


Rh 2 = WH2NH2, 

Rd2 - Wd2Nd2^ 

Rhd - WhdNhd- 

From, Fig. 2, distinct features are observed for various isotopes of 7 / 2 . Solid lines 
represents the results obtained from Monte Carlo method whereas dotted lines are for 
the results of Rate equation method. Difference between these two methods could be 
understood from Fig. 2. Monte Carlo method gives exact production rates due to the 
consideration of randomness. Here, we are computing the efficiency by considering the 
average production during the last few seconds. We find an interesting feature in the 
results of Monte Carlo method (Fig. 2abc). We have seen a dip in the efficiency profile 
of HD at the location where efficiency profile of D 2 having a peak. This feature is 
prominent in case of the Olivine grain (Fig. 2a). A decreasing feature in the efficiency 
profile of HD is visible at T = 11.5 Tf (Fig. 2a) whereas at the same location maximum 
efficiency is obtained for D 2 . This means that due to efficient production of D 2 , some 
D atoms are used up, which in turn decrease the production efficiency of HD. In case 
of Fig. 2b and 2c this feature is not well pronounced due to the lower production of 
£>2. 

In all the cases. Rate equation method shows lower efficiency (under production) 
for the production of H 2 but higher efficiency (over production) in case of HD and 
D 2 . In this context, we can compare Fig. 1 and Fig. 2b since in both the cases, we 
use amorphous carbon grain. Fig. 2b, shows that the production efficiency of HD and 
£>2 are higher in case of Rate equation method, which justifies the over production of 
HD and D 2 in Rate equation method (Fig. I). Since Monte Carlo method is more 
accurate, we recommend that this method must be used for all surface reactions. In 
case of Olivine grains, efficiency window for H 2 is wider and extended from 8 Tf to 
10 K. For HD, this window falls in between 9 - \A K. For D 2 , window is very narrow 
and peak is obtained at around 11.5 K. For amorphous Carbon grain (Fig. 2b), this 
window belongs to 11 K - 16 K, 15 K - 19 77, 21 K - 22 K respectively for H 2 , D 2 
and HD. From, Fig. 2c, efficiency window comes out to be at 11 - 16 K, 14 - 19 K 
and 14 - 17 K respectively for H 2 , D 2 and HD. 

So far, we discuss gas phase abundances of various isotopes of hydrogen molecules 
by considering thermal desorption from interstellar grains only. However, there should 
be another probability of desorption by which gas phase would be populated. This is 
called spontaneous desorption process, where some fraction of species could be evap¬ 
orated just after the formation. This happens due to the energy liberated during some 
reactions. To account for this feature in our model, we consider that a factor (1 - jj) 
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Figure 3: Difference between cases when spontaneous desorption factor is included (dashed line) and not 
included (solid line). 
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of Nh 2 , Nd 2 Nhd is lost to the gas phase. According to Katz et al. (1999), for an 
Olivine grain, = 0.33 may be used and for amorphous carbon grain, ^h 2 - 0.413 
could be used. Due to unavailability of experimental data, here, we use same fi values 
for HD and D 2 ifiH 2 - Here too, random numbers are generated for each 

newly formed H 2 , D 2 and HD and a fraction of species are allowed to populate the 
gas phase upon its production. In Fig. 3, we compare the efficiency window of H 2 , D 2 
and HD by considering the with (dashed line) and without spontaneous desorption term 
(solid line). Here we consider the Olivine grain for hh — 10^ cm“^ and - 0.1. While 
spontaneous desorptions are allowed, molecules are started to produce efficiently in gas 
phase at somewhat lower temperatures. For example, production efficiency of H 2 attain 
a peak value at 8 K when spontaneous desorption term is not considered (solid line of 
Fig. 3). With spontaneous desorption factor (dashed line of Fig. 3) this comes out to 
be at 7.5 K. For HD, solid line (no desorption) of Fig. 3 attain its maximum efficiency 
at 8.5 K while dashed line (considering spontaneous desorption) also attain a peak at 
8.5 K. But, while the spontaneous desorption term taken into account, HD started to 
produce at somewhat lower temperature. For example, a moderate efficiency at 7.5 K 
is achieved for the production of HD molecule while spontaneous desorption term is 
considered but at 7.5 K, efficiency for the production of HD is ~ 0 when spontaneous 
desorption term is not considered. The reason behind is that at low temperatures, ther¬ 
mal desorption time scales are much longer and spontaneous desorption factor allow 
some species to release from the grain surface and populate the gas phase. As a re¬ 
sults, surface species starts to populate the gas phase little bit lower temperatures while 
considering the spontaneous desorption factor. 

Following Chaabouni et al. (2012), sticking coefficient of H and D are calculated 
by using Eqn. 8. For a Silicate grain, Chaabouni et al. (2012) considered the sticking 
parameter (Sq) = 1 for the both H and D al Tq = 25 K and 50 K. From Eqn. 6, it is 
clear that the parameter is heavily dependent on temperature. In Eig. 4, a comparison 
between the consideration of unity sticking coefficient (normal case and used for all 
the cases) and the consideration by Chaabouni et al. (2012) is shown for an Olivine 
grain with tin = 10^ cm“^ and ro - 0.1. It is interesting to see that in the regime of 
our simulation, when we are considering sticking coefficients from Chaabouni et al. 
(2012), productions of HD and D 2 are insignificant. Only significant amount of H 2 is 
produced in this case. In between 5 - 13 K, sticking coefficient of H as calculated from 
Chaabouni et al., (2012) varies in the range of 0.95 - 0.80. In case of D, it varies in the 
range of 0.79 - 0.57. 

In Eig. 5(a-c), efficiency windows are shown when number of grain sites is varied. 
Here we consider an Olivine grain kept inside n/j = 10^ cm“^ and = 0.1. Eig. 5abc 
shows efficiency window for Olivine grain having (a) 2.5 x 10^ (Eig. 5a), (b) 4x10“* 
(Fig. 5b) and (c) 1.6 x 10^ (Fig. 5c) number of sites respectively. It is interesting to note 
that for larger grains, all molecules are producing efficiently. For smaller grains size 
(with accretion limit grain, having ~ 2.5 x 10^ sites), D 2 is not producing. However, 
for the larger grains (1.6 x 10^) it is producing. From Fig. 5(a-c), it is clear that as we 
are increasing grain size, efficiency for the formation of D 2 increases. 

It is discussed earlier that Monte Carlo method would provide a good estimation 
of recombination efficiency and Rate equation method often over or under estimate 
production rate. To overcome these discrepancies, Chakrabarti et al. (2006ab) defined 
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Figure 4: Efficiency window of a Olivine grain by considering sticking coefficient 1 (solid lines) and sticking 
coefficient from Chaabouni et al. (2012) (dashed line). 
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Ax = Ox/S, where Ax is effective recombination rate of a surface species to recombine 
with another surface species and is hopping rate as defined earlier. This is due to the 
fact that in two dimensions, a species (random walker) needs some numbers of steps 
which is linearly proportional to the number of distinct sites on the grain (Montroll & 
Weiss, 1965). If the grain having total S number of sites then it should be proportional 
to 5. In reality, S should be replaced by S “(f), where a(t) depends on various physical 
and chemical properties of grain (Chakrabarti et ah, 2006ab, Das et ah, 2008). In Rate 
equation, this linearity constant (a(t)) is taken to be one but in actual case this may 
not be the case. Rate equation method is very popular because it is economical to run 
in computer, whereas Monte Carlo method takes a long computational time even to 
handle a small chemical network. Chakrabarti et al. (2006ab) first proposed the idea 
for calculating various values of a(t) to modify Rate equation in such a way that that 
one could have an educated estimation of real scenario by using simple Rate equation 
method as well. Chakrabarti et al. (2006ab) studied formation of H 2 on interstellar 
grains and proposed values of a{t) under various physical circumstances. Das et al. 
(2008) studied the formation of Water and Methanol around the dense cloud following 
the same approach. In both the cases, it was assumed that a steady state could be 
reached at f —> 00 and at steady state a(f) = a. Here, we consider diffused cloud 
condition, where H and D are randomly accreting and producing H 2 , D 2 and HD. Due 
to inherent nature of the diffuse cloud condition, steady state condition may never be 
achieved (or achieved during very late stage of evolution). So a steady state value of a 
could be misleading. Here, to avoid any discrepancies, we calculate a factor and named 
it scaling factor {S /) defined to be the following: 




Number of surface species X at time t by Monte Carlo method 
Number of surface species X at time t by Rate equation method 


This factor defines number of surface species as predicted by Monte Carlo method with 
respect that obtained by Rate equation method. So after using Rate equation method, 
if we multiply the outcome by a scaling factor {S /) defined above, we may have a 
results as predicted by Monte Carlo method. In Fig. 6, we show variation of 5/ for 
H 2 , £>2 and HD for an Olivine grain kept at 9 K, n^ - lO'* cm^^, S =4x10"^ with 
number density of the cloud. For better illustration purpose, S is multiplied by 10 
and S is multiplied by 5. It is interesting to see that as we increase number density, 
5/ is increasing. Physical significance is that as we are going to higher density regime, 
grain surfaces are more and more populated and the surface species could easily find 
its reactant partner to react, so the production is enhanced. 

For low density case, production could be delayed due to unavailability of suitable 
reactant partner. Around high density region, production efficiency increases as well 
as 5/ increases. In Table 2, we present S / for all species {H, D, H 2 , D 2 and HD) 
for various sets of binding energies at temperatures where their production efficiency 
is maximum. For Amorphous Carbon grain, we provide results for T - \5 K and for 
intermediate energy values (set 3 energy values), T - W K. Similar trends could be 
observed for all sets of energies. 

In order to see the effects of the physical parameters on the formation of H 2 , HD 
and D 2 , we define another parameter /? and named it catalytic capability {jj) by follow¬ 
ing Chakrabarti et al. (2006ab) and Das et al. (2008), which measures efficiency of 
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sites. 


- HjiMC) 

- D2(MC) 

- HD(MC) 


Table 2: Values of 5/ for a wide parameter space. 


Type of grain 

Temperature (K) 

Hydrogen number density (cm 

^Jh 

SfHj 

^ fo 


^fHD 



50 

1.88 

1.20 

0.25 

0.029 

0.20 

Olivine 

9 

100 

1.84 

1.25 

0.24 

0.04 

0.20 



500 

1.86 

1.37 

0.245 

0.04 

0.20 



1000 

2.35 

1.40 

0.30 

0.04 

0.20 



5000 

0.84 

1.45 

0.06 

0.06 

0.20 



10000 

1.21 

1.45 

0.09 

0.06 

0.21 



50 

1.78 

1.41 

0.22 

0.03 

0.19 

Amoiphous Carbon 

15 

100 

1.79 

1.41 
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Figure 6: Variation of 5/ of H 2 , Do and HD with number density of the cloud for an Olivine grain kept at 
9 K. For the better visibility of the figure, S and S are multiplied by 10 and 5 respectively. 
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Figure 7: Variation of (5 with the accretion rate per site (lower label of X axis) or number density of the cloud 
(upper label of the X axis). 
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formation of H 2 , HD and D 2 on grain surface for a given pair of species residing on it. 
If 6 Nd 2 be the number of D 2 formed in St time then average rate of creation of D 2 per 
pair of deuterium atom is given by, 




1 SNp, 
2Nd St 


(15) 


We identify inverse of this rate as average formation rate and it is given by, 


Tf(t) = S>'^ 2 /Ad. 


Thus, 


Sf^‘>2(S = Ao/ < Aoiit) > 

This yields as a function of time: 

^02(1) = logiAol < Aoi(f) >)llog(S) 
Similarly for H 2 and HD, 

pH.it) = logiAn! < AhxH) »llog(S), 


Pnoif) - logiAnol < Anoiit) >)/log(S), 


(16) 

(17) 

(18) 

(19) 

( 20 ) 


where. 


and 


{Am(t)} = 


1 6Nh2 
2Nh St 


{AnDiit)} 


1 


SNh 


(Nh+Nd) St 

P(t) is also a time dependent parameter. At f —> 00 , /3(t) —> p. As time evolves, grains 
are populated by the species. Since, surface coverage is increasing, production should 
be faster due to decrease of reaction zone. But there should also be a blocking effect 
(Das et al. 2008), due to which surface species could be locked for a while, which in 
turn could delay production process. So, value of P(t) depends on surface coverage 
as well. After some transient time steps, we are having a steady state value of P(t). 
Here, we are considering at f —> 00 , j3(t) —> p. In Fig. 7, we show the value of p with 
respect to the accretion rate per site, p values are calculated during the last few steps 
of our simulations (~ 10*^ sec). As expected, with increasing accretion rate, surface 
coverage of species increases. In this situation, surface species needs to travel lesser 
number of steps to find one suitable reactant partner. Low value of p thus signifies 
faster production whereas higher value represents slower production rate. 

In Fig. 8, number of various surface species are shown with variation of r^. Obser¬ 
vational evidences suggest that elemental atomic DjH ratio (r/j) in an ISM would be 
~ 1.5 X 10“® (Linsky et al., 1995). Here, we consider an Olivine grain having 4 x 10“* 
sites kept at T - 9 K, hh - 10^ cm“^ and vary ro from 10“^ to 0.1. Due to random¬ 
ness of Monte Carlo method, here, we plot average numbers (average taken from last 
few steps after ~ 10* sec). For lower value of ro, production of D 2 is insignificant. 
Abundance of H and H 2 remain almost constant throughout simulation range of ro. 
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Figure 8: Valuation of the number of various surface species with the variation ro- 
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4. Conclusion 


In this paper, we mainly focused on production of H 2 , D 2 and HD on grain surfaces. 
Following are highlights of our results: 

• Production of these species highly dependent on the type of grains (i.e., interac¬ 
tion energies). We carried out our simulations using three types of binding energies. 

For Olivine grain, efficiency window is in between ^ K - 14 K, for amorphous carbon 
grain it is in between 11 K - 22 K and for binding energies as in LBH (set 3 energy 
values), this window is shifted to 11 /f - 19 /f. 

• We define a parameter S /, which solely depend on various physical and chemical 
properties of interstellar grains. Rate equation method often over or under estimates the 
production efficiency. To obtain more accurate production of these simple yet the most 
abundant molecules around ISM, this correction term should be considered in rectify¬ 
ing results from Rate equation method. For the sake of wider usage of our parameters, 
we provided a Table (Table 2) with various values of 5/ for a range of physical param¬ 
eters. 

• We computed another quantity /?, named catalytic capacity as used in earlier 
papers. This parameter shows a decreasing trend with increase in accretion rate. This 
implies production of surface species more and more favourable for high accretion 
regime. 

• Despite of low elemental abundances of atomic deuterium, several complex species 
are found to be heavily fractionated. Here, we vary initial DjH ratio to find out deu¬ 
terium fractionation of simple yet the most abundant species, namely, 7 / 2 . If we con¬ 
sider elemental D/7/ ratio of 10“^ (Linsky et ak, 1995) for a diffused cloud where all 
77 are in atomic form, production of D 2 found to be insignificant. 
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